library(tidyverse)
library(plotly)
library(gsignal)4. Modulazione di ampiezza
Esercizi sulla modulazione di ampiezza delle misure mediante analisi dello spettro e filtraggio offline.
1 Librerie
Quanto segue richiede questi pacchetti:
Se necessario installarli con install.packages() o mediante la GUI.
2 Funzioni di utilità
Per comodità definiamo alcune funzioni di utilità.
Generazione di un segnale sinusoidale con diverse armoniche, come definite in una tabella (pars) con le colonne w (ampiezza), f (frequenza) e phi (fase) per le varie frequenze (una riga per ogni armonica):
signal <- function(t, pars, rad = FALSE) {
stopifnot(is.data.frame(pars))
with(pars, {
if (!rad) {
phi <- phi/180*pi
f <- 2*pi*f
}
map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*sin(t*f[i] + phi[i] ))))
})
}3 Simulazione del segnale
Costruiamo l’uscita di uno strumento di misura, la sua modulazione in ampiezza, il rumore interferente e grafichiamo gli spettri.
plotly
Nel corso di statistica viene spiegato l’uso di GGPlot2. Quest’ultimo è studiato per l’efficacia e chiarezza di presentazione, per le quali è superiore a plotly.
È comunque possibile rendere i grafici GGPlot2 interattivi, passandoli alla funzione ggplotly().
Si noti tuttavia che avere molti grafici interattivi (e con molti punti) rallenta parecchio il caricamento della pagina nel browser: meglio evitarlo per i grafici “intermedi” e riservarlo ai grafici finali.
Definiamo alcuni parametri generici e la tabella delle armoniche:
# Parametri
fc <- 2000 # Frequenza di campionamento (Hz)
fm <- 10 # Frequenza del segnale di misura (Hz)
fr <- 15 # Frequenza del rumore (Hz)
fp <- 300 # Frequenza della portante (Hz)
m <- 0.8 # Indice di modulazione (0 < m < 1)
duration <- 10 # Durata del segnale (secondi)
# tabella delle armoniche del segnale di misura
pars_m <- tibble(
w = c(1, 0.5, 0.3),
f = c(fm, 2*fm, 4*fm),
phi = c(0, 0, 0)
)
pars_m %>% knitr::kable()| w | f | phi |
|---|---|---|
| 1.0 | 10 | 0 |
| 0.5 | 20 | 0 |
| 0.3 | 40 | 0 |
Generiamo il segnale di misura, perturbato da un rumore normale:
ym <- tibble(
t = seq(0, duration, by = 1/fc),
y = signal(t, pars_m),
yn = y + rnorm(length(t), 0, pars_m$w[1]/10)
) %>%
mutate(
f = 0:(length(t)-1)/max(t),
fft = fft(y)
)
ym %>%
mutate(intensity = Mod(fft) / length(t)*2,) %>%
ggplot(aes(x=f, y=intensity)) +
geom_line()
Generiamo il rumore interferente con il suo spettro. Parametri delle armoniche corrispondenti al rumore interferente:
pars_r <- tibble(
w = c(1, 0.5, 0.8),
f = c(fr, 3*fr, 9*fr),
phi = c(0, 40, 90)
)
pars_r %>% knitr::kable()| w | f | phi |
|---|---|---|
| 1.0 | 15 | 0 |
| 0.5 | 45 | 40 |
| 0.8 | 135 | 90 |
Aggiungo il rumore interferente alla tabella ym:
ym <- ym %>%
mutate(
y_in = signal(t, pars_r),
fft_in = fft(y_in)
)Ottengo il grafico dello spettro:
ym %>%
select(f, fft, fft_in) %>%
pivot_longer(
contains("fft"),
names_to = "signal",
values_to = "fft",
names_transform = ~ ifelse(.=="fft", "signal", "intf. noise")
) %>%
mutate(
intensity = Mod(fft) / length(t)*2
) %>%
ggplot(aes(x=f, y=intensity, color=signal)) +
geom_line()
Cerchiamo più possibile di utilizzare dati tidy, cioè una variabile per colonna, un’osservazione per riga. In questo modo realizzare grafici con serie multiple è particolarmente rapido (basta identificare la variabile di classificazione, cioè il nome della colonna che contiene la chiave di raggruppamento).
Per questo motivio, qui e di seguito useremo pivot_longer.
Aggiungo anche la modulazione:
# Aggiungo una mappa da nomi delle colonne -> nomi in legenda
fft_names=c(
"fft" = "signal",
"fft_in" = "intf. noise",
"fft_m" = "modulation"
)
# Aggiungo la modulazione
ym <- ym %>%
mutate(
y_m = (1 + m * yn) * cos(2 * pi * fp * t),
fft_m = fft(y_m)
)# Grafico dello spettro
(ym %>%
select(f, contains("fft")) %>%
pivot_longer(
-f,
names_to = "signal",
values_to = "fft",
names_transform = ~ fft_names[.]
) %>%
mutate(
intensity = Mod(fft) / length(t)*2
) %>%
ggplot(aes(x=f, y=intensity, color=signal)) +
geom_line() +
xlim(c(0, 500))) %>%
ggplotly()pivot_longer
Ridotto all’osso, l’ultimo grafico è (evito ggplotly per brevità):
ym %>%
select(f, contains("fft")) %>%
pivot_longer(-f) %>%
mutate(
intensity = Mod(value) / length(t)*2
) %>%
ggplot(aes(x=f, y=intensity, color=name)) +
geom_line()
Se non voglio usare pivot_longer() devo aggiungere tre colonne, una alla volta (stando attento a ripetere la formula correttamente: se scrivessi fft_m = Mod(fft_n) / length(t) * 2 otterrei un errore molto difficile da individuare!), con il modulo della FFT e nel plot devo aggiungere un layer alla volta (stando asttento a specificare correttamente il nome della serie come variabile color):
ym %>%
select(f, contains("fft")) %>%
mutate(
fft = Mod(fft) / length(t) * 2,
fft_in = Mod(fft_in) / length(t) * 2,
fft_m = Mod(fft_m) / length(t) * 2
) %>%
ggplot(aes(x=f)) +
geom_line(aes(y=fft, color="fft")) +
geom_line(aes(y=fft_in, color="fft_in")) +
geom_line(aes(y=fft_m, color="fft_m"))
4 Effetto interferente sul segnale modulato
Adesso che abbiamo il segnale modulato separato in frequenza dal rumore possiamo ‘esporlo’ al rumore interferente.
Osserviamo prima la modulazione:
ym %>%
ggplot(aes(x=t, y=y_m)) +
geom_line()
Combiniamo segnale e modulazione:
ym <- ym %>%
mutate(
y_m_in = y_m + y_in,
fft_m_in = fft(y_m_in)
)Ora ym è così composta:
ym %>% head() %>% knitr::kable()| t | y | yn | f | fft | y_in | fft_in | y_m | fft_m | y_m_in | fft_m_in |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.0000 | 0.0000000 | 0.1418529 | 0.0 | 0.00e+00+0.0000000i | 1.1213938 | 1.121394+0.0000000i | 1.1134823 | 3.829540+0.000000i | 2.2348761 | 4.950934+0.000000i |
| 0.0005 | 0.1004060 | 0.1040175 | 0.1 | 2.10e-06-0.0132478i | 1.1483848 | 1.121397-0.0073390i | 0.6366972 | 4.564120+5.526609i | 1.7850820 | 5.685517+5.519270i |
| 0.0010 | 0.2000641 | -0.0226319 | 0.2 | 8.30e-06-0.0265019i | 1.0386500 | 1.121407-0.0146797i | -0.3034221 | -3.960835-5.672208i | 0.7352280 | -2.839429-5.686887i |
| 0.0015 | 0.2982363 | 0.3224292 | 0.3 | 1.87e-05-0.0397689i | 0.8266720 | 1.121422-0.0220241i | -1.1963752 | 4.411363-12.118149i | -0.3697032 | 5.532785-12.140174i |
| 0.0020 | 0.3942043 | 0.4672451 | 0.4 | 3.33e-05-0.0530552i | 0.5637101 | 1.121445-0.0293739i | -1.1114244 | -5.362537+4.257905i | -0.5477143 | -4.241092+4.228531i |
| 0.0025 | 0.4872785 | 0.5019278 | 0.5 | 5.21e-05-0.0663672i | 0.3085893 | 1.121474-0.0367308i | 0.0000000 | 4.785527+5.336739i | 0.3085893 | 5.907001+5.300008i |
ym %>%
select(t, y_m, y_m_in) %>%
pivot_longer(-t, names_to = "series", values_to = "value") %>%
ggplot(aes(x=t, y=value, color=series)) +
geom_line() +
facet_wrap(~series, nrow=2) +
xlim(c(0, 0.5))
Osserviamo lo spettro:
fft_names["fft_m_in"] <- "mod. + intf."
ym %>%
select(f, contains("fft")) %>%
pivot_longer(
-f,
names_to = "signal",
values_to = "fft",
names_transform = ~ fft_names[.]
) %>%
mutate(
intensity = Mod(fft) / length(t)*2
) %>%
ggplot(aes(x=f, y=intensity, color=signal)) +
geom_line() +
xlim(c(0,500))
Estrarre il segnale in uscita originale a partire da quello affetto da rumore bianco, modulato ed affetto da rumore interfente ovvero da ym$y_m_in
4.1 In un colpo solo
Nella parte precedente, la tabella ym è stata costruita gradualmente, allo scopo di illustrare la composizione del segnale simulato. Guardando le operazioni nel loro complesso, la procedura è in realtà più semplice:
tibble(
t = seq(0, duration, by = 1/fc),
y = signal(t, pars_m) +
rnorm(length(t), 0, pars_m$w[1]/10), # segnale
fft = fft(y), # e sua FFT
y_in = signal(t, pars_r), # rumore interferente
fft_in = fft(y_in), # e sua FFT
y_m = (1 + m * y) * cos(2 * pi * fp * t), # modulante
fft_m = fft(y_m), # e sua FFT
y_m_in = y_m + y_in, # modulante + interferente
fft_m_in = fft(y_m_in) # e sua FFT
) %>%
# aggiungo la frequenza
mutate(
f = 0:(length(t)-1)/max(t),
) %>%
# prendo solo le colonne frequenza e fft
select(f, contains("fft")) %>%
# riorganizzo la tabella
pivot_longer(
-f,
names_to = "signal",
values_to = "fft",
names_transform = ~ fft_names[.]
) %>%
# aggiungo l'intensità
mutate(
intensity = Mod(fft) / length(t)*2
) %>%
ggplot(aes(x=f, y=intensity, color=signal)) +
geom_line() +
xlim(c(0,500))
5 OMETTERE E FAR FARE COME ESERCIZIO
Passo 1: Filtraggio e Demodulazione
Estraiamo il segnale modulato mediante un filtro passa-alto. Di conseguenza è rimasto solo il segnale modulato che va quindi demodulato per ottenere una stima del segnale in uscita originario
cutoff <- 2 * max(pars_r$f) / fc # Frequenza di taglio (2 volte la massima frequenza del rumore)
bf <- butter(4, cutoff, type = "high")
ym <- ym %>%
mutate(
y_m_f = filtfilt(bf, y_m_in) # filtfilt non introduce ritardo
)Passo 2: Demodulazione
Moltiplicazione per la portante:
ym <- ym %>%
mutate(
demod_raw = 2 * y_m_f * cos(2 * pi * fp * t) - 1
)Passo 3: Filtraggio Passa-Basso con filtro Butterworth
cutoff <- 2 * max(pars_m$f) / fc # Frequenza di taglio (2 volte la massima frequenza del segnale informativo)
bf <- butter(4, cutoff, type = "low")
ym <- ym %>%
mutate(
demod = filtfilt(bf, demod_raw)
)Confronto dei segnali:
(ym %>%
mutate(y_noisy = y + y_in) %>%
select(t, y, y_noisy, demod) %>%
pivot_longer(-t) %>%
ggplot(aes(x=t, y=value, color=name)) +
geom_line() +
coord_cartesian(xlim=c(0, 0.5))) %>%
ggplotly()Si noti come il segnale originario fosse affetto da rumore bianco, poi modulato e quindi sovrapposto ad un effetto interferente a tre armoniche. Alla fine del processo abbiamo ottenuto il segnale molto simile al segnale originario con una riduzione vistosa anche del rumore interferente